{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kriging in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "import geostat as geo\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Input data: the measurements\n",
    "x1 = np.array( [61,63,64,68,71,73,75] )        # Locations of the points (x1=x, x2=y)\n",
    "x2 = np.array( [139,140,129,128,140,141,128] )  \n",
    "v = np.array( [477,96,227,646,606,791,783] )   # Value of the measurements \n",
    "\n",
    "# The data in matrix form = one line per data point\n",
    "x = np.transpose( (x1,x2) )\n",
    "\n",
    "# The locations where we want to estimate \n",
    "x10 = np.array([68, 66, 69]); \n",
    "x20 = np.array([135, 132, 134]);\n",
    "xu = np.transpose( (x10,x20) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Input data: the measurements\n",
    "x1 = np.array( [61,63,64] )    # Loations of the points (x1=x, x2=y)\n",
    "x2 = np.array( [139,140,129] )  \n",
    "v = np.array( [477,96,227] )   # Value of the measurements \n",
    "\n",
    "# The data in matrix form = one line per data point\n",
    "x = np.transpose( (x1,x2) )\n",
    "\n",
    "# The locations where we want to estimate \n",
    "x10 = np.array([62, 63]); \n",
    "x20 = np.array([135, 133]);\n",
    "xu = np.transpose( (x10,x20) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot the data for visual check\n",
    "\n",
    "plt.plot(x1,x2,'+');   # Position des points\n",
    "plt.plot(x10,x20,'or') # Point où l'on veut estimer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ordinary kriging equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Defines the type of covariance function and parameters\n",
    "covar = geo.covariance(rang=10, sill=20,typ=\"spherical\")\n",
    "\n",
    "# Defines the size of the problem and agglomerate the data\n",
    "n_data = len(v)\n",
    "n_unknown = len(x10)\n",
    "n_total = n_unknown + n_data\n",
    "xall = np.concatenate((xu,x))\n",
    "\n",
    "# Computes the covariance for all pairs of nodes\n",
    "d = squareform( pdist(xall) )\n",
    "c = covar(d)\n",
    "\n",
    "# Assemble the Kriging matrix \n",
    "a = np.ones( (n_data+1,n_data+1) )\n",
    "a[0:n_data,0:n_data] = c[n_unknown:n_total,n_unknown:n_total]  \n",
    "a[n_data,n_data] = 0\n",
    "\n",
    "# Assemble the right hand side \n",
    "b = np.ones( (n_data+1,n_unknown) )\n",
    "b[0:n_data,0:n_unknown] = c[n_unknown:n_total,0:n_unknown]\n",
    "\n",
    "# Solve the Kriging system\n",
    "l = np.dot(np.linalg.inv(a), b)\n",
    "\n",
    "# Get the Kriging weight for each data point \n",
    "w = l[0:n_data,0:n_unknown]\n",
    "\n",
    "# Prepare the calculation of Kriging values \n",
    "v.shape=(n_data,1)\n",
    "\n",
    "# Computes the Estimated values by multiplying the kriging weights with values\n",
    "vo = np.dot( np.transpose(w), v)\n",
    "\n",
    "# Get the Lagrange parameter for each unknown point\n",
    "mu = l[-1,:] \n",
    "mu.shape=(n_unknown,1)\n",
    "\n",
    "# Computes the kriging variance for each unknown point\n",
    "so = np.diag( np.dot( np.transpose(w), b[0:n_data,:]) )\n",
    "so.shape = (n_unknown,1)\n",
    "so = so + covar._sill - mu\n",
    "\n",
    "print('Estimation:',vo)\n",
    "\n",
    "print('Variance:',so)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
